Fast followers?/Speed contagion: assessing the impact of Montreal F1 Grand Prix on high-speed ticketing rates (2000-2022)
Step 1. Weather variables. We defined the dates by year and availability (whether the event took place or not), days of the week, time windows, and pre and post-event spans, which we then linked to nearby meteorological stations. Also we linked with Collisions data.
Author
Andrés González Santa Cruz
Published
July 7, 2025
Code
# remove objects and memoryrm(list=ls());gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 835773 44.7 1627308 87 1176511 62.9
Vcells 1680897 12.9 8388608 64 3434063 26.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Adjuntando el paquete: 'dagitty'
The following object is masked from 'package:rio':
convert
Adjuntando el paquete: 'ggdag'
The following object is masked from 'package:stats':
filter
Adjuntando el paquete: 'kableExtra'
The following object is masked from 'package:dplyr':
group_rows
Cargando paquete requerido: Matrix
Adjuntando el paquete: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Adjuntando el paquete: 'lme4'
The following object is masked from 'package:rio':
factorize
This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
Cargando paquete requerido: carData
Adjuntando el paquete: 'car'
The following object is masked from 'package:purrr':
some
The following object is masked from 'package:dplyr':
recode
Cargando paquete requerido: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Adjuntando el paquete: 'brms'
The following object is masked from 'package:glmmTMB':
lognormal
The following object is masked from 'package:lme4':
ngrps
The following object is masked from 'package:stats':
ar
This is bayesplot version 1.12.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Adjuntando el paquete: 'bayesplot'
The following object is masked from 'package:brms':
rhat
This is loo version 2.8.0
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
- Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).
##
## Synth Package: Implements Synthetic Control Methods.
## See https://web.stanford.edu/~jhain/synthpage.html for additional information.
Warning: package 'weathercan' was built under R version 4.4.3
As of v0.7.2, the `normals` column in `stations()` reflects whether or not there
are *any* normals available (not just the most recent).
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
Registered S3 method overwritten by 'lfe':
method from
nobs.felm broom
Adjuntando el paquete: 'plm'
The following objects are masked from 'package:dplyr':
between, lag, lead
Adjuntando el paquete: 'PanelMatch'
The following object is masked from 'package:tidyr':
extract
The following object is masked from 'package:stats':
weights
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
Cargando paquete requerido: bsts
Cargando paquete requerido: BoomSpikeSlab
Cargando paquete requerido: Boom
Adjuntando el paquete: 'Boom'
The following objects are masked from 'package:brms':
ddirichlet, rdirichlet
The following object is masked from 'package:stats':
rWishart
Adjuntando el paquete: 'BoomSpikeSlab'
The following object is masked from 'package:stats':
knots
Cargando paquete requerido: zoo
Adjuntando el paquete: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Cargando paquete requerido: xts
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Adjuntando el paquete: 'xts'
The following objects are masked from 'package:dplyr':
first, last
Adjuntando el paquete: 'bsts'
The following object is masked from 'package:BoomSpikeSlab':
SuggestBurn
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
Adjuntando el paquete: 'forecast'
The following object is masked from 'package:brms':
ma
Code
#special repository indicated or the packageif(!require(weathercan)){install.packages("weathercan", repos =c("https://ropensci.r-universe.dev", "https://cloud.r-project.org")); library(weathercan) }#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_if(!require(bpmn)){devtools::install_github("bergant/bpmn")}
Warning in readLines(paste0(gsub("f1/", "", getwd()), "/key.txt")): incomplete
final line found on 'H:/Mi unidad/PERSONAL
ANDRES/UCH_salud_publica/pasantia/f1/key.txt'
API key has been set for the current session.
Collapse
We collapsed data by aggregating it at the day–year level for treatment and control areas. For each cell, we summed the number of high-speed tickets and documented vehicles and license holders, and averaged the minimum temperature, maximum temperature, while the 2-day-lagged precipitation and total daily precipitation were summarized by the median. We restricted post-outcome period to 7 days.
Code
#total_precip_median_lin lag_2_prec_median_lin_impif(identical(collisions_weather_corr_rect$lag_2_prec_median_lin_imp, collisions_weather_corr_rect$total_precip_median_lin)){warning("Identical lagged with not lagged variable")}collisions_weather_corr_rect_synth <- collisions_weather_corr_rect |>filter(date_num <=unclass(race_date) +7) |>group_by(tr_contr_sens, year, yday_corr, date) |>summarise(sum_velocidad =sum(velocidad, na.rm =TRUE),exp_veh =sum(vehicles_use_type, na.rm =TRUE), # Nº de vehículos expuestosexp_lic =sum(license_holders_sex_age, na.rm =TRUE), # Nº de conductoresmedian_lag_2_prec_median_imp =median(lag_2_prec_median_lin_imp, na.rm =TRUE),mean_min_temp_mean_lin =mean(min_temp_mean_lin, na.rm =TRUE),mean_max_temp_mean_lin =mean(max_temp_mean_lin, na.rm =TRUE),median_total_precip_median_lin =median(total_precip_median_lin, na.rm =TRUE),.groups ="drop" ) |>mutate( # ── aux variablesrate_veh = sum_velocidad / exp_veh *1e6, # not to modelrate_lic = sum_velocidad / exp_lic *1e6, # not to modeloff_veh =log(exp_veh), # offsetoff_lic =log(exp_lic) ) |>left_join(races, by =c("date"="race_date")) |>mutate(race_day =if_else(!is.na(year.y), 1L, 0L),id =as.numeric(factor(paste0(tr_contr_sens, year.x))) ) |>select(-year.y)if(identical(collisions_weather_corr_rect_synth$median_total_precip_median_lin, collisions_weather_corr_rect_synth$median_lag_2_prec_median_imp)){warning("Identical lagged with not lagged variable")}cat("Day of the race \n")collisions_weather_corr_rect_synth |>filter(race_day==1) |>reframe(min=yday_corr, max=yday_corr) |>print(n=40)############################################################## 1. “Only-treated, one series” ─────────────────────────## (collisions_weather_corr_rect_synth_tr)############################################################cat("Generate a summary of counts of every year in treated only\n")collisions_weather_corr_rect_synth_tr <- collisions_weather_corr_rect_synth %>%# start from the long tablefilter(tr_contr_sens =="treatment") %>%# keep treated rows onlygroup_by(tr_contr_sens, yday_corr) %>%# collapse across MRC × yearssummarise(sum_velocidad =sum(sum_velocidad, na.rm =TRUE),exp_veh =sum(exp_veh, na.rm =TRUE),exp_lic =sum(exp_lic, na.rm =TRUE),mean_median_lag_2_prec_median_imp =mean(median_lag_2_prec_median_imp, na.rm =TRUE),mean_min_temp_mean_lin =mean(mean_min_temp_mean_lin, na.rm =TRUE),mean_max_temp_mean_lin =mean(mean_max_temp_mean_lin, na.rm =TRUE),mean_median_total_precip_median_lin =mean(median_total_precip_median_lin, na.rm =TRUE),.groups ="drop" ) %>%mutate(rate_veh = sum_velocidad / exp_veh *1e6,rate_lic = sum_velocidad / exp_lic *1e6,off_veh =log(exp_veh),off_lic =log(exp_lic),year.x =2000, # dummy year / date just to keep structuredate =as.Date("2019-05-19"),race_day =0L,id =40# single ID for the synthetic treated series )############################################################## 2. “One treated - many controls” ─────────────────────## (collisions_weather_corr_rect_synth_one_tr_many_cntrls)############################################################# controls left at their original resolution (one ID per MRC–year)controls_long <- collisions_weather_corr_rect_synth %>%filter(tr_contr_sens !="treatment") %>%rename(mean_median_lag_2_prec_median_imp = median_lag_2_prec_median_imp,mean_median_total_precip_median_lin = median_total_precip_median_lin )collisions_weather_corr_rect_synth_one_tr_many_cntrls <-bind_rows( controls_long, collisions_weather_corr_rect_synth_tr # the single treated series from step 1)############################################################## 3. “Only-controls, one series” ───────────────────────## (collisions_weather_corr_rect_synth_cntr)############################################################collisions_weather_corr_rect_synth_cntr <- collisions_weather_corr_rect_synth %>%filter(tr_contr_sens !="treatment") %>%# keep controlsgroup_by(tr_contr_sens, yday_corr) %>%# collapse across all control IDssummarise(sum_velocidad =sum(sum_velocidad, na.rm =TRUE),exp_veh =sum(exp_veh, na.rm =TRUE),exp_lic =sum(exp_lic, na.rm =TRUE),# distributional summaries to keep a sense of variabilityp25_lag_2_prec_median_imp =quantile(median_lag_2_prec_median_imp, .25, na.rm =TRUE),p75_lag_2_prec_median_imp =quantile(median_lag_2_prec_median_imp, .75, na.rm =TRUE),mean_median_lag_2_prec_median_imp =mean(median_lag_2_prec_median_imp, na.rm =TRUE),p25_min_temp_mean_lin =quantile(mean_min_temp_mean_lin, .25, na.rm =TRUE),p75_min_temp_mean_lin =quantile(mean_min_temp_mean_lin, .75, na.rm =TRUE),mean_min_temp_mean_lin =mean(mean_min_temp_mean_lin, na.rm =TRUE),p25_max_temp_mean_lin =quantile(mean_max_temp_mean_lin, .25, na.rm =TRUE),p75_max_temp_mean_lin =quantile(mean_max_temp_mean_lin, .75, na.rm =TRUE),mean_max_temp_mean_lin =mean(mean_max_temp_mean_lin, na.rm =TRUE),p25_total_precip_median_lin =quantile(median_total_precip_median_lin, .25, na.rm =TRUE),p75_total_precip_median_lin =quantile(median_total_precip_median_lin, .75, na.rm =TRUE),mean_median_total_precip_median_lin =mean(median_total_precip_median_lin, na.rm =TRUE),.groups ="drop" ) %>%mutate(rate_veh = sum_velocidad / exp_veh *1e6,rate_lic = sum_velocidad / exp_lic *1e6,off_veh =log(exp_veh),off_lic =log(exp_lic),year.x =2000,date =as.Date("2019-05-19"),race_day =0L,id =1# single ID for the synthetic control ) %>%# rename the median variables so their names match the treated tablerename(mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin = mean_median_total_precip_median_lin )############################################################## Quick sanity checks############################################################# glimpse(collisions_weather_corr_rect_synth_tr)# glimpse(collisions_weather_corr_rect_synth_one_tr_many_cntrls)# glimpse(collisions_weather_corr_rect_synth_cntr)cat("Check imbalances\n")collisions_weather_corr_rect_synth_one_tr_many_cntrls |>group_by(id) |>reframe(max=max(yday_corr), min=min(yday_corr), n=n()) |>tail()invisible("Consolidate one treated and one control")collisions_weather_corr_rect_synth_cntr_tr <-rbind.data.frame(collisions_weather_corr_rect_synth_cntr,collisions_weather_corr_rect_synth_tr |># rename(collisions_weather_corr_rect_synth_tr, "median_lag_2_prec_mean_imp"= "mean_lag_2_prec_mean_imp", "median_total_precip_mean_lin"= "mean_total_precip_mean_lin")|> mutate(p25_lag_2_prec_median_imp=0, p75_lag_2_prec_median_imp=0, p25_min_temp_mean_lin=0, p75_min_temp_mean_lin=0, p25_max_temp_mean_lin=0, p75_max_temp_mean_lin=0, p25_total_precip_median_lin=0, p75_total_precip_median_lin=0))|>data.frame()
collisions_weather_corr_rect_quebec_synth <- collisions_weather_corr_rect |>filter(date_num <=unclass(race_date) +7) |>group_by(tr_contr, year, yday_corr, date) |>summarise(sum_velocidad =sum(velocidad, na.rm =TRUE),exp_veh =sum(vehicles_use_type, na.rm =TRUE), # Nº de vehículos expuestosexp_lic =sum(license_holders_sex_age, na.rm =TRUE), # Nº de conductoresmedian_lag_2_prec_median_imp =median(lag_2_prec_median_lin_imp, na.rm =TRUE),mean_min_temp_mean_lin =mean(min_temp_mean_lin, na.rm =TRUE),mean_max_temp_mean_lin =mean(max_temp_mean_lin, na.rm =TRUE),median_total_precip_median_lin =median(total_precip_median_lin, na.rm =TRUE),.groups ="drop" ) |>mutate( # ── aux variablesrate_veh = sum_velocidad / exp_veh *1e6, # not to modelrate_lic = sum_velocidad / exp_lic *1e6, # not to modeloff_veh =log(exp_veh), # offsetoff_lic =log(exp_lic) ) |>left_join(races, by =c("date"="race_date")) |>mutate(race_day =if_else(!is.na(year.y), 1L, 0L),id =as.numeric(factor(paste0(tr_contr, year.x))) ) |>select(-year.y)collisions_weather_corr_rect_quebec_synth_tr <- collisions_weather_corr_rect_quebec_synth %>%# start from the long tablefilter(tr_contr =="treatment") %>%# keep treated rows onlygroup_by(tr_contr, yday_corr) %>%# collapse across MRC × yearssummarise(sum_velocidad =sum(sum_velocidad, na.rm =TRUE),exp_veh =sum(exp_veh, na.rm =TRUE),exp_lic =sum(exp_lic, na.rm =TRUE),mean_median_lag_2_prec_median_imp =mean(median_lag_2_prec_median_imp, na.rm =TRUE),mean_min_temp_mean_lin =mean(mean_min_temp_mean_lin, na.rm =TRUE),mean_max_temp_mean_lin =mean(mean_max_temp_mean_lin, na.rm =TRUE),mean_median_total_precip_median_lin =mean(median_total_precip_median_lin, na.rm =TRUE),.groups ="drop" ) %>%mutate(rate_veh = sum_velocidad / exp_veh *1e6,rate_lic = sum_velocidad / exp_lic *1e6,off_veh =log(exp_veh),off_lic =log(exp_lic),year.x =2000, # dummy year / date just to keep structuredate =as.Date("2019-05-19"),race_day =0L,id =40# single ID for the synthetic treated series )controls_quebec_long <- collisions_weather_corr_rect_quebec_synth %>%filter(tr_contr !="treatment") %>%rename(mean_median_lag_2_prec_median_imp = median_lag_2_prec_median_imp,mean_median_total_precip_median_lin = median_total_precip_median_lin )collisions_weather_corr_rect_quebec_synth_one_tr_many_cntrls <-bind_rows( controls_quebec_long, collisions_weather_corr_rect_quebec_synth_tr # the single treated series from step 1)collisions_weather_corr_rect_quebec_synth_cntr <- collisions_weather_corr_rect_quebec_synth %>%filter(tr_contr !="treatment") %>%# keep controlsgroup_by(tr_contr, yday_corr) %>%# collapse across all control IDssummarise(sum_velocidad =sum(sum_velocidad, na.rm =TRUE),exp_veh =sum(exp_veh, na.rm =TRUE),exp_lic =sum(exp_lic, na.rm =TRUE),# distributional summaries to keep a sense of variabilityp25_lag_2_prec_median_imp =quantile(median_lag_2_prec_median_imp, .25, na.rm =TRUE),p75_lag_2_prec_median_imp =quantile(median_lag_2_prec_median_imp, .75, na.rm =TRUE),mean_median_lag_2_prec_median_imp =mean(median_lag_2_prec_median_imp, na.rm =TRUE),p25_min_temp_mean_lin =quantile(mean_min_temp_mean_lin, .25, na.rm =TRUE),p75_min_temp_mean_lin =quantile(mean_min_temp_mean_lin, .75, na.rm =TRUE),mean_min_temp_mean_lin =mean(mean_min_temp_mean_lin, na.rm =TRUE),p25_max_temp_mean_lin =quantile(mean_max_temp_mean_lin, .25, na.rm =TRUE),p75_max_temp_mean_lin =quantile(mean_max_temp_mean_lin, .75, na.rm =TRUE),mean_max_temp_mean_lin =mean(mean_max_temp_mean_lin, na.rm =TRUE),p25_total_precip_median_lin =quantile(median_total_precip_median_lin, .25, na.rm =TRUE),p75_total_precip_median_lin =quantile(median_total_precip_median_lin, .75, na.rm =TRUE),mean_median_total_precip_median_lin =mean(median_total_precip_median_lin, na.rm =TRUE),.groups ="drop" ) %>%mutate(rate_veh = sum_velocidad / exp_veh *1e6,rate_lic = sum_velocidad / exp_lic *1e6,off_veh =log(exp_veh),off_lic =log(exp_lic),year.x =2000,date =as.Date("2019-05-19"),race_day =0L,id =1# single ID for the synthetic control ) %>%# rename the median variables so their names match the treated tablerename(mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin = mean_median_total_precip_median_lin )collisions_weather_corr_rect_quebec_synth_cntr_tr <-rbind.data.frame(collisions_weather_corr_rect_quebec_synth_cntr,collisions_weather_corr_rect_quebec_synth_tr |>mutate(p25_lag_2_prec_median_imp=0, p75_lag_2_prec_median_imp=0, p25_min_temp_mean_lin=0, p75_min_temp_mean_lin=0, p25_max_temp_mean_lin=0, p75_max_temp_mean_lin=0, p25_total_precip_median_lin=0, p75_total_precip_median_lin=0))|>data.frame()
Test series structure
Code
acf(collisions_weather_corr_rect_synth_tr$rate_veh*10, lag.max =50, main ="Treatments\nACF:high-speed tickets per 1MM vehicles")pacf(collisions_weather_corr_rect_synth_tr$rate_veh*10, lag.max =50, main ="Treatments\nPACF:high-speed tickets per 1MM vehicles")
ACF plots, treated
ACF plots, treated
We noticed a partial autocorrelation of 10 days in the treated.
Code
acf(collisions_weather_corr_rect_synth_cntr$rate_veh*10, lag.max =50, main ="Controls,\nACF:high-speed tickets per 1MM vehicles")pacf(collisions_weather_corr_rect_synth_cntr$rate_veh*10, lag.max =50, main ="Controls,\nPACF:high-speed tickets per 1MM vehicles")
ACF plots, controls
ACF plots, controls
We noticed a partial autocorrelation of 25 days in the controls.
Code
x_vals <-time(forecast::msts( collisions_weather_corr_rect_synth_cntr$rate_veh*10,start =c(2019, 1),seasonal.periods =c(7)))forecast::autoplot(forecast::mstl(forecast::msts(collisions_weather_corr_rect_synth_tr$rate_veh*10,start =c(2019, 1),seasonal.periods =c(7)), robust =TRUE))+ggtitle("Treated")+scale_x_continuous(breaks =seq(2019, 2025, by =1/7)[seq(1, 43, by =7)], # where the ticks sitlabels =seq(-35, 7)[seq(1, 43, by =7)] # what the ticks say ) +geom_vline(xintercept =seq(2019, 2025, by =1/7)[which.min(abs(seq(-35, 7)))], # rel_day 0linetype ="dashed", colour ="black")forecast::autoplot(forecast::mstl(forecast::msts(collisions_weather_corr_rect_synth_cntr$rate_veh*10,start =c(2019, 1),seasonal.periods =c(7)), robust =TRUE))+ggtitle("Controls")+scale_x_continuous(breaks =seq(2019, 2025, by =1/7)[seq(1, 43, by =7)], # where the ticks sitlabels =seq(-35, 7)[seq(1, 43, by =7)] # what the ticks say ) +geom_vline(xintercept =seq(2019, 2025, by =1/7)[which.min(abs(seq(-35, 7)))], # rel_day 0linetype ="dashed", colour ="black")
Weekly decomposition of series
Weekly decomposition of series
The decomposition shows that differences between treated and control groups emerge in the trend, not in weekly seasonality, with treated show a downward trend, while controls show an upward trend. Any modelling or graphical comparison should therefore adjust for the weekly effect and focus on the diverging long-term trajectory.
We visualize series differentiating by treated or control cities.
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = sum_velocidad)) +geom_line(aes(group = id, colour =ifelse(id==40,"Treated", "Controls"))) +scale_colour_brewer(palette ="Dark2", name ="MRC") +geom_vline(xintercept =30, linetype ="dashed", color ="black", linewidth=1) +labs(x ="Days (-30 to 5 days post-treatment)", y ="Number of High‑speed Tickets (sum)") +theme_minimal() +theme(panel.grid.minor =element_blank(), legend.position ="bottom")
Daily High‑speed tickets per Treated/Control cities
Code
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = (sum_velocidad/exp_veh)*1e6)) +geom_line(aes(group = id, colour =ifelse(id==40,"Treated", "Controls"))) +scale_colour_brewer(palette ="Dark2", name ="MRC") +geom_vline(xintercept =30, linetype ="dashed", color ="black", linewidth=1) +labs(x ="Days (-30 to 5 days post-treatment)", y ="Number of Daily high-speed collisions per 1MM vehicle counts (sum)") +theme_minimal() +theme(panel.grid.minor =element_blank(), legend.position ="bottom")
Rate of daily high-speed collisions per 1MM vehicle counts by Treated/Control status
Code
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_min_temp_mean_lin)) +geom_line(aes(group = id, colour =ifelse(id==40,"Treated", "Controls"))) +scale_colour_brewer(palette ="Dark2", name ="MRC") +geom_vline(xintercept =30, linetype ="dashed", color ="black", linewidth=1) +labs(x ="Days (-30 to 5 days post-treatment)", y ="Min temperature") +geom_ribbon(data = collisions_weather_corr_rect_synth_cntr,aes(x = yday_corr, y = mean_min_temp_mean_lin,ymin = p25_min_temp_mean_lin,ymax = p75_min_temp_mean_lin),inherit.aes =FALSE,fill ="grey70", # o un color que combine con tu paletaalpha =0.3 )+theme_minimal() +theme(panel.grid.minor =element_blank(), legend.position ="bottom")
Min temperature by Treated/Control status
Code
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_max_temp_mean_lin)) +geom_line(aes(group = id, colour =ifelse(id==40,"Treated", "Controls"))) +scale_colour_brewer(palette ="Dark2", name ="MRC") +geom_vline(xintercept =30, linetype ="dashed", color ="black", linewidth=1) +labs(x ="Days (-30 to 5 days post-treatment)", y ="Max temperature") +geom_ribbon(data =subset(collisions_weather_corr_rect_synth_cntr_tr, id !=40),aes(x = yday_corr, y = mean_max_temp_mean_lin,ymin = p25_max_temp_mean_lin,ymax = p75_max_temp_mean_lin),inherit.aes =FALSE,fill ="grey70", # o un color que combine con tu paletaalpha =0.3 )+theme_minimal() +theme(panel.grid.minor =element_blank(), legend.position ="bottom")
Max temperature by Treated/Control status
Code
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_median_total_precip_median_lin)) +geom_line(aes(group = id, colour =ifelse(id==40,"Treated", "Controls"))) +scale_colour_brewer(palette ="Dark2", name ="MRC") +geom_vline(xintercept =30, linetype ="dashed", color ="black", linewidth=1) +labs(x ="Days (-30 to 5 days post-treatment)", y ="Precipitations \n(annual mean of city-level medians, where each city median\nis computed from the medians of its weather stations)") +geom_ribbon(data =subset(collisions_weather_corr_rect_synth_cntr_tr, id !=40),aes(x = yday_corr, y = mean_median_total_precip_median_lin,ymin = p25_total_precip_median_lin,ymax = p75_total_precip_median_lin),inherit.aes =FALSE,fill ="grey70", # o un color que combine con tu paletaalpha =0.3 )+theme_minimal() +theme(panel.grid.minor =element_blank(), legend.position ="bottom")
Total precipitations by Treated/Control status
Code
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_median_lag_2_prec_median_imp)) +geom_line(aes(group = id, colour =ifelse(id==40,"Treated", "Controls"))) +scale_colour_brewer(palette ="Dark2", name ="MRC") +geom_vline(xintercept =30, linetype ="dashed", color ="black", linewidth=1) +labs(x ="Days (-30 to 5 days post-treatment)", y ="Precipitations (2-day lagged)\n(annual mean of city-level medians, where each city median\nis computed from the medians of its weather stations)") +geom_ribbon(data =subset(collisions_weather_corr_rect_synth_cntr_tr, id !=40),aes(x = yday_corr, y = mean_median_lag_2_prec_median_imp,ymin = p25_lag_2_prec_median_imp,ymax = p75_lag_2_prec_median_imp),inherit.aes =FALSE,fill ="grey70", # o un color que combine con tu paletaalpha =0.3 )+theme_minimal() +theme(panel.grid.minor =element_blank(), legend.position ="bottom")
Total precipitations (lagged 2) by Treated/Control status
Analysis plan
Code
library(DiagrammeR)gr <-grViz("digraph study_design {graph [layout = dot, rankdir = TB]node [shape = rectangle, style = filled, color = LightSkyBlue, fontname = Helvetica]start [label = 'Analysis Start']adelanto [label = 'Advance two days (D)']no_adelanto [label = 'No advance\n(from race day) (D_off)']quebec_only [label = 'Quebec only as control (D)']quebec_sherbrooke [label = 'Quebec and Sherbrooke\nas control (D_sens)']exp_3days [label = 'Exposure 3 days\nafter the race (D)']exp_7days [label = 'Exposure 7 days\nafter the race (D7)']# edgesstart -> {adelanto no_adelanto}adelanto -> {quebec_only quebec_sherbrooke}quebec_only -> {exp_3days exp_7days}}", width =1000,height =1400)gr
Code
unlink(paste0(getwd(),"/_figs/analysisplan_files"), recursive =TRUE)htmlwidgets::saveWidget(gr, paste0(getwd(),"/_figs/analysisplan.html"))webshot::webshot(paste0(getwd(),"/_figs/analysisplan.html"), paste0(getwd(),"/_figs/analysisplan.png"),vwidth =300, vheight =300*1.2, zoom=10, expand=100) # Test w/ diff coords:top, left, width, and height
LaZerte SE, Albers S. weathercan: Download and format weather data from environment and climate change canada. The Journal of Open Source Software. 2018;3(22):571. https://joss.theoj.org/papers/10.21105/joss.00571.